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Abstract. We present a high order parameter-robust numerical method for a sys¬ 
tem of {M > 2) coupled singularly perturbed parabolic reaction-diffusion problems. 
A small perturbation parameter e is multiplied with the second order spatial deriva¬ 
tives in all the equations. The parabolic boundary layer appears in the solution of 
the problem when the perturbation parameter e tends to zero. To obtain a high or¬ 
der approximation to the solution of this problem, we propose a numerical method 
that employs the Crank-Nicolson method on an uniform mesh in time direction, 
together with a hybrid finite difference scheme on a generalized Shishkin mesh in 
spatial direction. We prove that the resulting method is parameter-robust or e- 
uniform of second order in time and almost fourth order in spatial variable, if the 
discretization parameters satisfy a non-restrictive relation. Numerical experiments 
are presented to validate the theoretical results and also indicate that the relation 
between the discretization parameters is not necessary in practice. 


Key words. Singular perturbation, Parabolic reaction-diffusion problems. Cou¬ 
pled systems. High order compact scheme, Crank-Nicolson method. Parameter ro¬ 
bust method. Generalized Shishkin mesh. 


1. Introduction 


We consider the following system of (M > 2) coupled singularly 
perturbed parabolic react ion-diffusion problems 
3u 

(1) L^u ■= -^ + Lx,eU = f, (x, t) E D :=nx (0, T] = (0,1) x (0, T], 

(2) u{0,t) = 0, u{l,t) = 0, Vf G [0,T], 

(3) rt(x,0) = 0, Wx E fl. 

The spatial differential operator is dehned by 

0 \ / aii(a:) 


/ —e— 


^x,e 


V 0 


—£ 


dx^ 


-1-A, with A = 


/ 






o-mm{x) 


where £ is a small parameter that satishes 0 < e -C 1. Denote the 
boundaries of the domain D by F := Fq IJ Ti) with Fq = {(x, 0)|x E D} 
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and Fi = {{x,t)\x = 0,1, t G [0,T]}. We assume that the coupling 
matrix A = {aij{x))MxM satishes the following positivity conditions at 
each X E Q 

(4) ttij ^ 0, i 7^ j, 

M 

(5) an >0, ^ ajj > /3* > 0, i = 1,..., M. 

1=1 

If (5) is not satished directly, we consider the transformation u{x,t) = 
u{x,t) exp(—/Sot) with (3o > 0 (sufficiently large) in order to transform 
the diagonal entries such that (5) holds. Also, we assume that sufficient 
regularity and compatibility conditions hold among the data of the 
problem (l)-(3) such that the exact solution u G In the 

analysis we assume the following compatibility conditions (see 0) 

The numerical analysis of singular perturbation problems has always 
suffered from serious difficulties due to the boundary layer behavior 
of the solution when the perturbation parameter becomes small. Re¬ 
cent years have witnessed substantial progress in the development of 
layer adapted meshes to design a special class of numerical methods, 
so called parameter-robust numerical methods, that converge uniformly 
with respect to the perturbation parameter (see m)- Parameter- 
robust numerical methods based on htted meshes, particularly the 
Shishkin meshes gained popularity because of their simplicity and ap¬ 
plicability to more complicated problems in higher dimensions, see 

(6) for more details. Several numerical studies for coupled system 
of singularly perturbed reaction-diffusion problems are considered in 
[iQ] , [H] , [l2] , [IS] , [IE] and the references therein. 

To solve the system of two coupled singularly perturbed parabolic 
reaction-diffusion problems with the distinct small perturbation pa¬ 
rameters in each equations, Gracia and Lisbona [5] proposed a uni¬ 
formly convergent numerical method by using the classical backward 
Euler scheme in time and the central difference scheme in spatial di¬ 
rection, and proved that the error bound is 0{At In^ N) with 

the assumption N~'^ ^ CAt, 0 < q < 1. High order numerical meth¬ 
ods have always been an interest for the numerical community as they 
provide good numerical approximations with low computational cost. 
Recently, Clavero et ah [1] gave an attempt to design a high order 
uniformly convergent numerical method for solving the system of two 
coupled singulary perturbed parabolic reaction-diffusion problems with 




PARABOLIC REACTION-DIFFUSION SYSTEMS 


3 


the distinct small perturbation parameters in each equations. To in¬ 
crease the order of uniform convergence, the authors in [1] considered 
the Crank-Nicolson method on an uniform mesh in time direction and 
central difference scheme on a standard Shishkin mesh in spatial direc¬ 
tion, and proved that the error bound is 0((At)^ -|- In^ N) with 

the assumption ^ CAt, 0 < q < 1. To our knowledge this is 
the only high order parameter-robust numerical method is available in 
the literature for solving parabolic react ion-diffusion system (l)-(3). In 
the present paper, our objective is to integrate the available techniques 
for high order approximations (eg. |1] and 0), to design a high or¬ 
der parameter-robust numerical method for solving parabolic reaction- 
diffusion system (l)-(3). For a high order approximation, we consider 
the Crank-Nicolson method on an uniform mesh in time, together with 
a hybrid scheme which is a suitable combination of the fourth order 
compact difference scheme and the standard central difference scheme 
on a generalized Shishkin mesh in spatial direction. It can be seen 
that the combination of Crank-Nicolson method in time direction with 
hybrid scheme in spatial direction does not satisfy the discrete max¬ 
imum principle except if the restrictive condition At ^ C{L/N)‘^ is 
imposed. In this article, we follow the approach of Clavero et ah [2] 
to overcome this difficulty. First, some auxiliary problems are con¬ 
sidered which permits to prove appropriate bounds for local error of 
the Crank-Nicolson method. Then the uniform convergence analysis 
of the scheme used to discretize these auxiliary problems is discussed. 
Finally, using the recursive arguments and the uniform stability of the 
totally discrete scheme, we claim that the present method is uniformly 
convergent of second order in time and almost fourth order in spatial 
variable. It should be noted here that in the theoretical proof we as¬ 
sume the totally discrete scheme operator satisfy the uniform stability 
as a conjecture in Section 5. As so far it is an open problem to prove the 
uniform stability of totally discrete scheme theoretically (see also ID)- 
While in the support we presented the numerical tables (Tables 2,4 
and 6) that shows the spectral radius of the totally discrete operator is 
strictly less than one, independent of e and discretization parameters, 
in Section 6. 

This paper is arranged as follows. In Section 2, a priori bounds on 
the solution of (l)-(3) and its derivative are constructed. The time 
semidiscretization using the Crank-Nicolson method and its local con¬ 
sistency error is given in Section 3. In this section we also discuss the 
asymptotic behavior of the solution of semidiscretized problems and 
their spatial derivatives. In Section 4, the generalized Shishkin mesh is 
given and the spatial semidiscretization with a hybrid scheme which is 
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a suitable combination of the fourth order compact difference scheme 
and the central difference scheme is described on generalized Shishkin 
mesh for the set of stationary singularly perturbed problems studied in 
Section 3. It is also proved that the spatial semidiscretization is almost 
fourth order uniformly convergent on generalized Shishkin mesh. In 
Section 5, semidiscretization steps are combined to give the total dis¬ 
cretization and its uniform convergence is proved. The numerical ex¬ 
periments are conducted to demonstrate the efficiency of the proposed 
method in Section 6. Finally, conclusions are included in Section 7. 

Notations; In the remaining parts of the paper, C and (7 = (7(1,..., 1)^ 
are the generic positive constant independent of £ and discretization 
parameters. Dehne n^mifuj^tCj, and |n| = 

(|ni|,..., \vm\Y'■ We consider the maximum norm and it is denoted 
by ||.||h, where is a closed and bounded set. For a real valued func¬ 
tion V G C{H) and for a vector valued function v = (ui,... ,vm)'^ G 
C{H)^, we dehne 

||n||H = max |n(a;)| and ||n||H = max{||ni|..., \\vm\\h}- 

xGH 

li H = Q, we drop H from the notation. The analogous discrete 

_ g 

maximum norm on the mesh is denoted by 11. | Us . For any function 

_ _ “N 

g G (7(r2), gi is used for g{xi)-, if g G C{Vt)^ then gj = g{xi) = 
{gi^i,... ,gM,i)'^- For simplicity, we use Latj, for L{Nq). If Nq = N, we 
drop N as subscript from the notation L^. 

2. Properties of the exact solution 

Following the technique of Theorem 1 in [5], we can show that the 
operator in (1) satishes the following maximum principle. 

Lemma 2.1. Let y G fl . Let y{x, 0) > 0 on and 

y(0, t) > 0, y{l, t) > 0 on [0, T]. Then L^y > 0 in D implies y^ 0 on 
D. 


An immediate consequence of Lemma 2.1 is the following stability 
result. 

Lemma 2.2. Let u be the solution of (l)-(3). Then 

To obtain the bounds on the solution u of (l)-(3), the variable x 
is transformed to the stretched variable x dehned hy x = xfy/e, this 
results that (l)-(3) transformed as 
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(6) 


^ ^ ^ 

L^u := — + = f, {x,t) e 


(7) 

where 


Lx,e 


/ _ 92 
9x2 


V 0 


u{x,t) = 0, {x,t) e Te, 

0 \ ^ / aii(^) ••• aiM(^) 

+A, with ^ = I : ■ •. : 

-^ / \ aMi{x) ... aMM{x) 


Ds = fie X (0,T] = (0, l/y/e) X (0,T] and Fg is its boundary analogous 
to r. Here the differential equation (6) is independent of e. Using the 
standard local estimate for the solution of system of time dependent 
partial differential equations (see |H]), we obtain the bounds on the 
solution of (6)-(7) and its derivative. On returning in term to the 
original variable x and using ||m|| 7 j ^ O, obtained from Lemma 2.2 
with e-uniform boundedness of /, yields the following result. 


Lemma 2.3. Let u he the solution of (l)-(3). Then it satisfies 

^ for0^i + 2j < 6. 

T) 

In the following result, we derive sharper bounds on the derivatives 
of u to show that the large values seen in Lemma 2.3 do in fact decay 
rapidly as one moves away from the boundary F. 


dx^dP 


Lemma 2.4. Let u be the solution of (l)-(3). Let (3 G (0,/3*) be arbi¬ 
trary but has a fixed value. Then there exists a constant C, independent 
of e, such that 
( 8 ) 

< c(l + £-”*/2(exp(-a:,//?/£) + exp(-(l - x)^&/e)3j 


d'^u{x, t) 


for {x, t) E D and m = 0,..., 6. 

Proof Fix (3 G (0,/3*) and set Pm{x) = 1 -|- e~"^P{exp{—x^/(3/e) -|- 
exp(—(1 — x)-\/pJe)). The proof is by mathematical induction. The 
bound ([H]) for m = 0 follows from Lemma 2.2. Assume that ([H]) holds 
for m = 0,..., 12 — 1, 1 ^ z/ ^ 6. We now prove ([S]) for m = v. Letting 
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note that 



y{x,0) = 0 


A 


{y—^) d^u 
dx^ 


in D, 

in f2, 
in (0,T], 


where boundary conditions follow from Lemma 2.3. From the induc¬ 
tive hypothesis, it is clear that ^ CP^_i{x). Applying the 

maximum principle with the barrier function CPy{x)^ we obtain the 
required result, i.e., for {x,t) G D 


d'^u{x, t) 
dx^ 


^ C {l + e ''/^(exp(-a;vW^) + exp(-(l - x) vW^))) • 


This proves the lemma. 


□ 


Now a special decomposition of the exact solution u into a regu¬ 
lar part V and a layer part w can be obtain as follow. Set x* = 


4 .Je/l 3 ln(l/Ve). 

Dehne for each j G 

{ 1 ,- 

.., n} and 

(x, t) E D 

(9) 

f 4 






E 

(x - x*y ^ 


for 0 ^ X 

^ x*,t E [ 0 ,T]; 

Vj{x, t) = < 

u =0 

Uj{x 

4 



for X* ^ a 

■ ^ 1 — X*, t G [ 0 , T] 


E 

^ u =0 

{x-x*y 

\t) 

for 1 — X* 

^ X ^ 1 , t G [ 0 ,T], 


and Wj{x, t) = Uj{x, t) — Vj{x, t). Then Lemma [2^ and the choice of x* 
yields, for s = 0,..., 6, (cf. Linss [S]) 


(10a) 


d^Vj{x, t) 

dx^ 


^ C{1 + 


(10b) 

^ ^exp(-x\//3/£) -1- exp(-(l - x) vW^)) • 

It should be noted here that this decomposition does not, in general, 
satisfy L^v = f and L^w = 0 . 


d^Wj{x, t) 

(9x* 


3. The time semidiscretization 

We introduce the time semidiscretization of (l)-(3) by using the clas¬ 
sical Crank-Nicolson method, with constant time step At on uniform 
mesh vj = {nAt, 0 ^ n ^ T/At}. The time semidiscretization is given 
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by 

( 'iP = u{x, 0) = 0 , 

(11) (j+f i^ju-+i = (j-fx,,,)w- + f(r+r+'), 

[ (0) = (1) = 0, for n = 0,1,..., T/At - 1, 

where is the approximation of the exact solution u of (l)-(3) at the 
time level = nAt, u = 0,1, 2,..., T/At and /"■ = f{x, tn). 

To study the consistency of (ITT]) , we dehne the following auxiliary 
problem 

( 12 ) 

r := (/ + fX.,.)«"+' = (/ - f X,,,)n(a;,t„) + f (/- +/-+i), 

\ ^”+^(0) = w"+^(l) = 0, 

where is the approximation to u{x, tn+i)- Let e„+i(a;) = u{x, tn+i) — 
u^^^{x) be the local truncation error of fllll) and it satishes the follow¬ 
ing lemma. 

Lemma 3.1. If 

^ C, {x,t) E D, 0 ^ i ^ 3, 

then the local error associated to the scheme m satisfies 

|e^_i_i(a;)| ^ C{Atfi, x E 

Proof. The results follows from the arguments given in [2]. □ 

Now we prove that the asymptotic behavior of the solution of the 
semidiscretize problem flT^ and its spatial derivative have essentially 
the same asymptotic behavior that the solution of a stationary system 
of coupled singularly perturbed react ion-diffusion problems. Using the 
approach Clavero et ah [1], such feature is given by the following lemma. 


d'^u{x, t) 
dt^ 


Lemma 3.2. Let be the solution of 17^) . Then it satisfies 
(13) 

^ C^^l + e~^^‘^{exp{-x^/^) -Fexp(-(1 - x)^/fije)^ , 

where 0 ^ /c ^ 6 and C is a constant independent of e and At. 

Proof. Let us hrst start by studying the behaviour of , that means 
the result flT^ for A; = 0. As the data f is e-uniformly bounded, 
\u{x,tn)\ ^ C and \Lx^eU{xfin)\ ^ C] similar to |T2], the operator 
(/ -|- ^La- e) satishes a maximum principle and using this it follows 
that 

^ c. 


d^u 

dx^ 
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To prove the result flT^ for the derivatives of , we introduce the 
following auxiliary function 

0”+^ = - u{x,tn)), 

which is the solution of the following boundary value problem 

r (/+f = -2L,,eu{^, tn) +r+ r^\ 

( 14 ) 

[ </)^+i(o) = 0, (jf^+^iX) = 0. 

Using \L^^euX,tn)\ = \f{x,tn) - ^(x,4)| ^ C and |0''+^(O)| ^ 
(7, ^ C with the maximum principle for (/ + 

get 


^ C. 


\r*'\ 

Next we write the problem flT^ as 

^)X"-L^,Mx,tn)+r+r^^, 


T -^n+l 

Lx,e'^ 


(15) 


u 


72+1 


( 0 ) = 0 , = 0 


From ^ C, it can be seen that the right side of flT^ is £- 

uniformly bounded. Using this with +”"'"^1 ^ C we get 


(16) 


^2-n+l 


'It 




^ Ce 


-1 


From (fT6|) and using the mean value theorem argument as used in 
we obtain 


(17) 


dx 


€ Ce-^l\ 


On differentiating flT^ with respect to x, we define 
1, 2, are the solutions of boundary value problems 

{(/ + fn,.)cr' = 9.(n, 

( 18 ) 1 

y cr‘(o) = si cr‘(i) = si 

where |s:^| ^ i = 1,2, j = 0,1 and using Lemma 2.4 

\gi{x)\ ^ C'(l + £“*/^(exp(-xvW^)+exp(-(l-x)vW^)))i * = 1,2- 
Now taking the barrier function as 
Xx) = Ciil + x) + C' 2 £“*''^(exp(-xvW^) + exp(-(l - x)yW^)) 


•n+1 


dx'^ ’ 


I = 


and for sufficiently large value of Ci and C 2 using the maximum prin¬ 
ciple for (J -|- ^Lx^e)i we deduce that flT^ is true for /c = 1, 2. 
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Now to prove the bound ffT^ for higher value of k, we follow similar 
arguments given in [T]. This proves the lemma. □ 


Next we dehne the Shishkin-type decomposition for the solution of 
semidiscretize problem flT^ . This type of decomposition has been dis¬ 
cussed earlier in Linss [5], for scalar singularly perturbed boundary 
value problem. To dehne this, let x* = {4:^/eJ^)\n{l/^/e). Similar 
to (9), for each x G hi and k = we set for 

X G [x*, 1 — X*] and extends to a smooth function dehned on Q 
and dehne ~ x G hi. Then the results of Lemma 

3.2 and the choice of x* implies the following decomposition of 
(cf. Linss [9]) 


Lemma 3.3. Let be the solution of / f7^) . Then it can be repre¬ 
sented as , where the regular part satisfies 


(19) 


dx"* 




and the layer part satisfies 


( 20 ) 


d'^w 


n+l 


dx^'' 


< 


for 0 ^ m ^ 6, k 
and At. 


Ce (^exp(-xvW^) + exp(-(l - x) vW^)) > 

= 1,..., M, and C is a constant independent of e 


The above lemma shows that the solution of flT^ is decom¬ 
posed into a sum of regular part and layer 

part , that is, it can be written as = 

This decomposition is said to be a Shishkin-type de¬ 
composition (not a standard Shishkin decomposition) as it does not in 
general satisfy (/-h = {I - ^L^,e)v{x,tn) + 

and {I = {I — ^Lx^i;)w{x,tn), as these additional prop¬ 

erties are not needed in the error analysis of present method. 


4. The spatial semidiscretization 


In this section, hrst, we construct a generalized Shishkin mesh S{L) 
to discretized the domain hi := [0,1] by using a suitable mesh generat¬ 
ing function % as described in im. Dehne the transition parameter 


a = min 


1 

4’ 


a^y/sL 




where cto(> is a positive constant and L = L(iV) the value of L 

with iV elements that satisfy In(lniV) < L ^InN and 



(21) 
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_ g 

The mesh points of generalized-Shishkin discretized domain are 
given by Xj = %{j/N), j = 0,1,, N/2, and by symmetry = 

1 — Xj, j = 0,1,..., N/2, where X G C'^p, 1/2] and dehned as 

f Aat, for t G [0,1/4]; 

(22) X(f) = I 

[ p{t - 1/4)^ + 4(T(f - 1/4) + a, for t G [1/4,1/2]. 
Here the coefficient p is determined by X(l/2) = 1/2. 

_ Q 

Note that the mesh 12jy is uniform in [0, a] and [1 — a, 1], and it 
changes smoothly in the transition points {a, 1 — a}. However, the 
mesh width hj = Xj+i — Xj, for j = N/i ,..., 3iV/4, satishes (see [T7] ) 

^ N-^ max X'(t) ^ 

[{i-l)/N,{i+l)/N] 

^ ’ h.+i - hj ^ iV-2 max X" (t) ^ CN-\ 

[(i-l)/N,(i+l)/N] 

We shall let hmax = max hj,j = 1,2,..., N, and by symmetry it is 
easy to verify that hmax = /iv /2 = /iA/ 2 +i- 

4.1. The hybrid scheme. We introduce a hybrid scheme to dis¬ 
cretize the set of stationary coupled system of singulary perturbed 
reaction-diffusion problem flT^ on the generalized Shishkin mesh 12^. 
The hybrid scheme is a combination of the fourth order compact dif¬ 
ference scheme (where the coefficients g/’s and r/’s of the scheme are 
determined so that the scheme is exact for the polynomials up to degree 
four and satisfy the normalization conditions = 1, i = 

1,2,..., N — 1, k = 1,2,..., M ) and the central difference scheme, 
and is given by 


(24) 

—n+1^ 

U ]. = 

[Q?\ for i = 1,2, 

...,iV-l, 

(26) 

-n+1 

Uo =0, 

U = 0, 


where 




(26) 

/ R{U^+^) 

+ ^Q{0‘12U2~^^) + • 

.. + f Q(aiMfo;^+') 

—-n+1^ 

|i„t/ ]i- 

i?(f/-+i) + fQ^a2iU^^^) + . 

.. + f Q(a2Mfo]^+') 


rNP) + 

fQ(aMifor') + ... 
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( 27 ) 


(28) 

with 



/ Qifi) \ 
Qif^) 

V Qif^) ) 


f R{Vk,i) 

1 QiVk,) 



Vk,i-1 + + 'f'i’~^Vk,i+l 

Vk,i-1 + Qi'^Vk,i + (li’~^^k,i+l 




7 (Xi) = u{Xi,tn) + ^{-[L^,eU]{Xi,tn) +f{Xi,tn) + f {Xi, tn+l)) ■ 


The coefficients ’ ,i = 1,... N — 1, k = 1,2,..., M, * = —, c, + are 
given by 
(29) 

fimArr, + + s)) 

M{„k,- 
2 

At ( -2e 


k,— 

r.’ = 

k.c 
i 


K,C At/ K,— I fc,C I \ fc,— k,+ I 1 

fc,^ 


-2e I I 2 


The coefficients gf’*) * = 1, • • •, -/V — 1, A; = 1, 2 ..., M, * = —, c, + are 
defined in two different ways. 


(t) For the mesh points located in (0, r) U (1 — r, 1), the coefficients 
g7*,t = {l,...,iV/4-l}U{3iV/4 + l,...,iV-l},A) = l,2,...,M, * = 
—, c, +, are given by 


(30) 





1 

T2‘ 


{ii) For the mesh points located in [r, 1 — r], depending on the re¬ 
lation between /tmax and e, the coefficients q^'*^ t = 1,..., iV — 1, k = 
1,2, ...,M, * = —,c,-|- are defined in two different cases. Define 
Ofcfc = Ofcfc 2 /At for fc = 1 , 2 ,..., M. 


In the first case, when T^^axl loo ^ where 7 is a positive con¬ 
stant independent oie and At, the coefficients gf’*, i = iV/4,..., 3A/4, k 
1 , 2 ,..., M, * = — , c, -f, are given by 

(31) 

6yhj 6 -1- /ij_|_i) 

While in the second case, when loo > £, where 7 is a 

positive constant independent of e and At, the coefficients gf’*,i = 
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A^/4,..., 3iV/4, /c = 1, 2,..., M, * = —,c,+, are given by 
(32) gf- = 0, = 1, gf + = 0. 

The above definition of coefficients gf’s and r^’s show that the scheme 
(24)-(25) is defined by the fourth order compact difference scheme 
within the boundary layer region (0, r) U (1 — r, 1). While in the regu¬ 
lar region [r, 1 — t], the scheme (24)-(25) is defined by a modified high 
order non-equidistant difference scheme when loo ^ ^ and 

is defined by the central difference scheme when qh^axl loo > 
This means that the scheme (24)-(25) considers the high-order approx¬ 
imation only when the local mesh width is small enough to give non¬ 
positive off-diagonal entries while at all other mesh points the central 
difference scheme is used. This combination leads to the following 
lemma. 

Lemma 4.1. Let 7 = 1/6 and No be the smallest positive integer such 
that 

max{4(To(||afcfc||oo + 2/Af)/3} < (No/Ln^)^, 

k 

where = L{Nq) defined in (21). Then, for any N > Nq, the 
discrete operator defined by (2f)-(25) is of positive type. 

Proof. Firstly, for Xi G (0, r) U (1 — r, 1), the fourth order com¬ 
pact difference scheme is considered in this region. The condition 
max{4cro(||afcfc||oo + 2/Af)/3} < {No/L^^fi for any N > Nq, where 

k 

Ltvo = L{No) as defined in (21), with the the coefficients , r^*, * = 
—, c, -|-, = 1, 2,..., M, defined by fl2^ - fl30|) and the assumption (4)- 

(5), concludes the lemma. 

Secondly, for Xi G [t, 1 — r] when T^^axl l®fcfc| |oo > the central dif¬ 
ference scheme is considered. Hence the proof is trivial. 

While in the opposite case, for Xi G [r, 1 — r] when 7h^axl l®fcfc| |cx) ^ 
we decompose [r, 1 — r] := [r, X 7 V/ 2 ] U [a:Ar/ 2 ,l — r] and study the 
sign of coefficients g/s in these two cases, separately. First, when 
Xi G [r, a:Ar/ 2 ] the coefficients gf’"'" is clearly non-negative on general¬ 
ized Shishkin mesh while the coefficient gf’~ will be non-negative 
when 2hi — /i^+i ^ 0 for A/4 ^ i ^ A/2. The assertion is triv¬ 
ially true for i = A/2 because of uniform mesh (at symmetry). For 
A/4 ^ i ^ A/2 - 1, we can write it 

hi ^ hi, 

follows if (cf. (23)) 

AX'((i - 1)/A) ^ X"((i + 1)/A), 
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that is, if 

w{z) = 3pz'^ — 6pz — 12p + 4criV^ ^ 0, 

where z = i — 1 — iV/4 ^ —1. It is not hard to verify that the dis¬ 
criminant of the quadratic function w is non-positive if 4criV^ ^ 15p. 
Since qh^axl loo ^ ^ dehnitely implies 4criV^ ^ 15p, it follows that 
w{z) >0 for all 2 ;. This completes the proof. 

Similar to this, we can prove ^ 0 for N/2 ^ i ^ 3iV/4. Thus 
the condition qh^axl loo ^ ^ with the coefficients gf’*, rf’*, * = 
— ,c,+, k = 1,2,, M, dehned by (12^ .131) and the assumption (4)- 
(5), concludes the lemma. 

□ 

Remark 4.2. It can he seen that the scheme (24)-(25) on standard 
Shishkin mesh does not satisfy the above Lemma f.l because the coeffi¬ 
cients gf’* are not always non-negative at the transition points, due to 
the fact that the standard Shishkin mesh is very anisotropic in nature. 
While if we consider the scheme (24)-(25) on the generalized Shishkin 
mesh 12^ then the coefficients gf’* are always non-negative. This is 
used in the proof of Lemma 4-L 


Using the Lemma 4.1, the discretization operator dehned by (24)- 
(25) is of positive type and it satishes the following discrete comparison 
principle. 

Lemma 4.3. (Discrete Comparison Principle) Let V and W be two 
mesh functions and satisfy V\i > ^ W\i, i = 1,2,..., N—1, Vq > 

Wq and Vn > W^, then Vi > W) i = 0,1,..., N. 


Using the above discrete comparison principle we obtain the follow¬ 
ing discrete stability estimate. 

Lemma 4.4. (Discrete Stability Estimate) Let V be the mesh function 
with Vq = Vn = 0. Then 

wn^cwil^n 

where C is independent of N, At and e. 

Let r-n-i-i(xi) be the truncation error associated to the scheme (24)- 
(25) and is dehned by 
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Lemma 4.5. Let be the solution of and U be the ap¬ 
proximate solution of the spatial discretized scheme (2f)-(25). Let the 
hypothesis of Lemma 4-1 be satisfied. Then the global error satisfies 

^ CA^(L/iV)^ 

with the assumption that L~^ < CAt, where C and C are a positive 
constants independent of N, At and e. 

Proof. If r = 1/4, then the mesh Djsi is uniform, that is, N~^ is very 
small respect to e and therefore a classical analysis can be used to prove 
the convergence of the scheme. So, in the analysis we only consider the 
case r = a^y/eL. 


The truncation error estimate r-n+i(a:j) of the scheme (24)-(25) on the 
generalized Shishkin mesh Hjv is discussed in the following cases. 


(/) When Xi G (0, r) U (1 — r, 1), we have hi = hj+i = Aaoy/eN ^L. 
By Taylor expansions we obtain 


(33) 


ir. 


n+1 


fa:,) I ^ CeAth'f 




dx^ 




Using hi = Aaoy/eN ^L and 


rjri+l 


dx^ 


^ Ce it follows that 


(34) I [r(4x)], I ^ CAt{L/N)\ 

{II) When Xi G [r, 1 — r ], according to the decomposition = 
, split the truncation error into two parts to obtain 


(35) I r2n+i(xi) I < I r^„+i(xi) I +1 r^n+i(xi) |. 

For the mesh points located in [r, 1 — t], depending on the relation be¬ 
tween hmax and e, the scheme (24)-(25) is dehned by the combination 
modihed high order non-equidistant difference scheme and the central 
difference scheme. The error analysis for both cases are given as follows. 


{i) For the case qh^axl|oo ^ £, suppose g G C®[0,1]^, then by 
Taylor expansions we obtain 

(36) |r 3 ,fc(a:j)| ^ CeAt{Pk,i + Qk,i + Rk,i), ior k = 1,..., M, 

where 

Pk,i = {hi+l—hi) Qk,i ~ I I fl'i ^ I I pi-i.Xi+i]; 

Using fl2^ and flT^ . we obtain the bound of the truncation error with 
respect to the regular part 
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(37) |r^n+i(xi) I ^ CAtiV-^ 

Again using (1^ and fl20|) . we obtain the bound of the truncation error 
with respect to the layer part 

(38) I rg^™+i(a:i) | ^ |5,| 

For Xi e [r, 1 - r], 

||_g III- , < p{—^N/A-l\fFli) g(“(I“*3JV/4-|-l)-\/^A) — 2e(“^-'V/4-l vWe) 

1111 [Xi_i \ ~ 

= ^ where r = doV^L, do > 4/v^. 

Then e~^ ^ L/N leads to 

(39) 

Using ([39]) in (|3H|) with 7^LxlU ^ we get 

(40) \r^n+^{xi)\^ CAt{L/N)\ 

On combining fl57|) and (140]) with dOO]) . we obtain 

(41) I T^n+i{xi) I < C'At(L/A)^, for 7^Lxl|ofcfc|U ^ £• 

(a) Now, for the case 7h^axl loo > £, suppose g E ^^[0,1]^, then 
by Taylor expansions we obtain 

(42) |rg^fc(a:i)| ^ CeAt{Yk^i + Z^,*), for fc = 1,..., M, 
where 

Using (1231) and (llOp . we obtain the bound of the truncation error with 
respect to the regular part 

(43) |r^„+i(a:i) I ^ CeAtA-^. 

Now using the condition 7h^axl |cx) > we obtain 

(44) |r~n+i(xi)| ^ CA-^ 

Note that in (44) the term At disappears from the bound for the 
error associated with the regular part; this fact is important in order 
to impose the relation between the discretization parameters At and A. 

To estimate the error with respect to the layer part , suppose 
g E C^[0,1]^, then using 

I rg(a:j) I ^ CeAt^g ||[3:i_i,a:i+i], 
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we have 

|rsi-+i(xi)| ^ 

Using fl3^ . we have 

(45) \V^^+i{xi)\^C^t{LlNf. 

Combining (14^ and (1451) in (]35|1 with the assumption such that < 
CAt, we obtain 


(46) I T^n+i{xi) I ^ CAt{L/NY, for 7^maxlU > £• 

On combining the case (/) and case (//), we obtain the truncation error 
estimate for the scheme (24)-(25) on the generalized Shishkin mesh Dj^ 
and it is given by 

(47) |r~.+i(xi) I ^ C'At(L/iV)^ 

Therefore, from the truncation error estimate fITT)) and the uniform 
stability result given in Lemma 4.4, we conclude the lemma. □ 

5. Total discretization 


In this section, we write the total discretization by combining the 
time semidiscretization and spatial semidiscretization to compute the 
approximate solution of (l)-(3) and after that we prove that the resul¬ 
tant scheme is uniformly convergent of second order in time and almost 
fourth order in spatial variable. Concretely, the numerical approximate 
of u{xi, nAt) for i = 1,..., A and n = 0,1,..., T/At, are obtained 
by the following totally discrete scheme 


r C/? = 0, = £„m(x.,0), i = 0(l)Af, 

(48) I [£7= |gF”|„ for * = l(l)Af - 1, 

[7(^+1 = 0, C/”+i = 0, forn = 0,l,...,T/At-l, 


where 

(49) 



/ i?(f/r+i) + fQ(ai2Ur') + - + yQ(aiMU"+') \ 

i?(t/-+i) + ^Q^a2iU^+^) + ... + fQ{a2MU2+^) 

\ + fQ{aMiU^^^) + -. + fQ{aMM-iU^t\) J ■ 


[QF^V ■■= 


At 


Ql-Ff) \ 
Q{F?) 

q(F2,) y 


f "(i.) = tJ” + —(-£7 E/” + /(ij, („) +/(!., f„+i)), 


(50) 






PARABOLIC REACTION-DIFFUSION SYSTEMS 


17 


for z = l(l)iV — 1, n = 0,1,..., T/Af — 1, 




Un+l 
2 - * 


u: 


At 


+ f{Xi,tn) +f{Xi,tn+l), 


and 

(51) 


f RiVk,) 
I Q{Vk,i) 



Vk,i-1 + A ’ + A ’"*'14,1-1-1, 

14,i-l + <ii ’ 14,i + 4 


The coefficients gf’*, * = —, c, + and * = —, c, +, /c = 1, 2,..., M 
are defined as in Section 4. 


Theorem 5.1. Let u be the exact solution of (l)-(3) and let {h/”"*"^} 
he the numerical solution of the scheme (fS). Under the hypothesis of 
Lemma f.l, the global error u{x, tn+i) — C/”"*"^ at the time tn+i satisfies 

l|,i(xi,4„+i) - < c((Atf + (i/iv)*). 

with the assumption that L~^ < CAt, where C is a positive constant 
independent of N, At and e. 

Proof. The global error u{xi,tn+i) — of the totally discrete 

scheme at the time tn+i can be split in the form 

(52) ^ 

u{XiUn+l) - 4 {u{XiUn+l) “ {Xi)) + {vP""^{Xi) - ) 

+(uT - c^r^)- 

On combining the result from the Lemma 3.1 and Lemma 4.5 with 

(52) , we obtain 

(53) 

||i,„ $C((A«)UA«(L/Af)‘‘) + ||t/, 

-n-l-1 _i_i 11 -''^+1 _i_i 

To bound the term 11 , we consider that U — U 

can be written as the solution of one step of (48) with starting value 
u{xi,tn) — Uf, taking the source term f equal to zero together with 
zero boundary conditions. Then it follows that 

(54) t/i - t/”+‘ = J?A,(«(i.,(„)- t/”), 

where is a linear operator, called the transition operator associated 
to the totally discrete scheme (48). Using this with (53) we obtain a 
recursive argument as 

n 

(55) ||«(Xj,i„+i)-t/“+‘ « Cj]||fl”-'-||i,„((A()nA((L//V)‘‘). 

k=l 
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To get the required result for the uniform convergence of totally discrete 
scheme a sufficient condition is that 

By assuming the uniform boundedness condition on power of discrete 
transition operator R]\f with (53) (see the Remark 5.2 below) and the 
hypothesis of Lemma 4.5 that L~‘^ ^ CAt, we conclude the main results 
of this section. □ 

Remark 5.2. We assume here the uniform boundedness condition as 
a conjuncture holds for the transition operator R^, as the theoretical 
proof of this is an open problem so far in the literature. Some partial re¬ 
sults in this direction can be obtained by using a result by Palencia [T4] . 
but for the present problem this would require an e-uniform estimate of 
the resolvent of the spatial operator L^^e- Here due to lack of available 
theoretical result in this direction we assume the uniform boundedness 
of the power of discrete transition operator Rj^ as a conjuncture. For 
the support of this conjuncture we show some numerical evidence for 
the spectral radius of Rjq. From the numerical results of the Tables 2, 
4, and 6 we observe that the spectral radius of R^ is strictly less than 
one and it stabilize as the singular perturbation parameter e becomes 
small. 

Remark 5.3. Theorem 5.1 proves almost fourth order uniform conver¬ 
gence of the method in spatial variable under the relation L~^ ^ CAt. 
Nevertheless, from the numerical point of view in Section 6, this con¬ 
dition is an artificial relation that we never needed in the experiments. 
Note that this relation appeared when we prove the convergence of the 
regular components in regular region, see eg. (44) Section 4- 

6. Numerical experiments 

The proposed method is implemented on three test examples. In 
all the cases we begin with total number of nodal points iV = 64 and 
the time step At = 0.5. The maximum error at the nodal points is 
calculated for the different values of e and N. 

Example 1: Consider the following system of two coupled singularly 
perturbed parabolic problem 

dUi d'^Ui \ 2/1 \2 

+ {2 + x)ui - {I + x)u 2 = x^{l-xY, 

dU2 d‘^U2 , r /1 \ 2/1 \2 

-^-£-^^ + (e + 1)U2 - {I + x)ui = x^{l-xY, 

for {x,t) e (0,1) X (0,1], with the initial-boundary conditions 

u{x, 0) = 0, x G (0,1), 
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m(o, t) = u{i, t) = 0, t e [0, 1], 
and the exact solution is not known. 


Example 2: Consider the following system of three coupled singularly 
perturbed parabolic problem 

+ 3ui - {1 - x)u 2 - {I - x)u 3 = 16x^{l-xf, 

dU2 d‘^U2 / . N „ q 

-2u,-U3 = t , 

^ “ ^“2 = 16a:^(l - 

for {x,t) G D, with the initial-boundary conditions 

u{x, 0) = 0, X G (0,1), 


^(0, t) = u{l, t) = 0, t G [0,1], 
and the exact solution is not known. 


Example 3: Consider the following system of two coupled singularly 
perturbed parabolic problem 

dui d‘^Ui 


dt 

dU2 


— e 


dx'^ 

d‘^U2 


+ 2ui -U2 = 1, 
2 u 2 - Ml = 1, 


dt dx"^ 

for (x,f) G (0,1) X (0,1], with the initial-boundary conditions 

m(x, 0) = 0, X G (0,1), 

M(0,t) = ^(1,^) = 0, t G [0,1], 
and the exact solution is not known. 


As the exact solution is not known for these examples, we estimate 

the maximum nodal error, Es^N,At = max e^’^*(i,nAt), for different 

’ ’ Vi,72 

values of e and N where = |Mjv(xj,f„) — u^{xi,tn)\- We 

use a variant of the double mesh principle, assume Uf^{xi,tn) denotes 
the numerical solution at the nodal point (xj, on the tensor product 

_5 

mesh of the generalized Shishkin mesh with iV -|- 1 nodal points 
in spatial direction and a uniform mesh of step size At in time direc¬ 
tion, and U]\f{xi,tn) denotes the numerical solution at the nodal point 
{xi,tn) on the tensor product mesh {(xj,t„)} that contains the mesh 
points of the original mesh and their midpoints. 
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Table 1 . Maximum error and numerical rate of conver¬ 
gence of the present method with uniform step size At in 
time direction and generalized Shishkin mesh S{L) with 
L = L* in spatial direction for the Example 1. 


£ = 2-^ 

N = 64 iV = 128 N = 256 
At = 0.5 At = 0.5/4 At = 0.5/4^ 

iV = 512 N = 1024 
At = 0 . 5/43 ^ q_5/44 

k=A 

8.82E-04 

4.42E-05 

2.71E-06 

1.69E-07 

1.06E-08 


4.32 

4.03 

4.00 

4.00 


8 

4.19E-04 

2.57E-05 

1.60E-06 

l.OOE-07 

6.26E-09 


4.03 

4.00 

4.00 

4.00 


12 

3.91E-04 

2.43E-05 

1.52E-06 

9.52E-08 

5.95E-09 


4.01 

3.99 

4.00 

4.00 


16 

3.88E-04 

2.41E-05 

1.52E-06 

9.51E-08 

6.00E-09 


4.01 

3.99 

4.00 

3.99 


20 

3.86E-04 

2.41E-05 

1.52E-06 

9.48E-08 

5.93E-09 


4.00 

3.99 

4.00 

4.00 


24 

3.86E-04 

2.41E-05 

1.52E-06 

9.48E-08 

5.93E-09 


4.00 

3.99 

4.00 

4.00 


28 

3.86E-04 

2.41E-05 

1.52E-06 

9.48E-08 

5.93E-09 


4.00 

3.99 

4.00 

4.00 


32 

3.86E-04 

2.41E-05 

1.52E-06 

9.48E-08 

5.93E-09 


4.00 

3.99 

4.00 

4.00 


EN,At 

8.82E-04 

4.42E-05 

2.71E-06 

1.69E-07 

1.06E-08 

pN 

4.32 

4.03 

4.00 

4.00 



In the standard way, we estimate the classical convergence rate, for 
each hxed e, by 


jsf 

Pe = - 


ln(-£'£,2A,At/4) 


In 2 


and the parameter-robust convergence rate by 
N ln(-E'jv,At) — ln(-E'2Ar Af/4) 

^ =-hr2-^ 


where EN At = max E^ n At. 

’ V£ ’ ’ 
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Table 2. Spectral radius of the transition operator 
for the Example 1. 


e = 2-^ N = U iV = 128 iV = 256 N = 512 N = 1024 
At = 0.5 At = 0.5/4 At = 0.5/4^ At = 0.5/4^ At = 0.5/4^ 


k=4: 

0.40350 

0.80792 

0.94825 

0.98681 

0.99670 

8 

0.57817 

0.87118 

0.96616 

0.99143 

0.99785 

12 

0.59223 

0.87965 

0.96849 

0.99203 

0.99800 

16 

0.59810 

0.88169 

0.96905 

0.99217 

0.99804 

20 

0.59953 

0.88219 

0.96919 

0.99221 

0.99805 

24 

0.59989 

0.88234 

0.96923 

0.99222 

0.99805 

28 

0.59997 

0.88235 

0.96923 

0.99222 

0.99805 

32 

0.59999 

0.88235 

0.96923 

0.99222 

0.99805 


Using L < hi N instead of IniV; this means we are trying to bring 
the point xi closer to x = 0 and this provides the higher density of 
the mesh points in the layers. The motivation for this is the fact that 
the better performance of the mesh S{L) can be governed by the high 
density of mesh points in the layers. The smallest value of L is chosen 
to be L* = L*{N) which satishes 

e-^* = L*/N. 

For the different values of N and £, Table 1, Table 2, and Table 3 
represent the maximum error and the classical rate of conver¬ 

gence of the present method for the Example 1, Example 2, and 
Example 3, respectively. The last two rows in each of the tables (Table 
1, Table 3, and Table 5) represent the maximum error with respect to 
each nodal point for all value of e, that is ] and the parameter- 

robust numerical rate of convergence . 

To show the numerical evidence for the uniform stability of the tran¬ 
sition operator R^, we calculate the spectral radius of i^Ar for different 
value of N, At and e. Table 2, Table 4, and Table 6 display the spectral 
radius of this operator for the Example 1, Example 2, and Example 3, 
respectively. We clearly observe that the spectral radius for all value 
of N, At and e is always strictly less than one. Moreover, we observe 
that the spectral radius stabilized for the small value of singular per¬ 
turbation parameter e. This stabilization of spectral radius for small 
value of e indicates the uniform stability of the operator R^. 
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Table 3. Maximum error and numerical rate of conver¬ 
gence of the present method with uniform step size At in 
time direction and generalized Shishkin mesh S{L) with 
L = L* in spatial direction for the Example 2. 


£ = 2-^ 

N = 64 iV = 128 N = 256 
At = 0.5 At = 0.5/4 At = 0.5/4^ 

iV = 512 N = 1024 
At = 0 . 5/43 ^ q_5/44 

k=A 

4.48E-02 

4.09E-03 

2.17E-04 

1.35E-05 

8.44E-07 


3.45 

4.23 

4.01 

4.00 


8 

4.44E-02 

3.58E-03 

1.97E-04 

1.22E-05 

7.63E-07 


3.63 

4.18 

4.01 

4.00 


12 

4.43E-02 

3.54E-03 

1.95E-04 

1.21E-05 

7.58E-07 


3.65 

4.18 

4.01 

4.00 


16 

4.43E-02 

3.54E-03 

1.95E-04 

1.21E-05 

7.58E-07 


3.65 

4.18 

4.01 

4.00 


20 

4.43E-02 

3.54E-03 

1.95E-04 

1.21E-05 

7.58E-07 


3.65 

4.18 

4.01 

4.00 


24 

4.43E-02 

3.54E-03 

1.95E-04 

1.21E-05 

7.58E-07 


3.65 

4.18 

4.01 

4.00 


28 

4.43E-02 

3.54E-03 

1.95E-04 

1.21E-05 

7.58E-07 


3.65 

4.18 

4.01 

4.00 


32 

4.43E-02 

3.54E-03 

1.95E-04 

1.21E-05 

7.58E-07 


3.65 

4.18 

4.01 

4.00 


EN,At 

4.48E-02 

4.09E-03 

2.17E-04 

1.35E-05 

8.44E-07 

pN 

3.45 

4.23 

4.01 

4.00 



Observe that the data in Example 3 does not satisfy the zeroth order 
compatibility conditions at the nodal points (0, 0) and (1, 0). Moreover, 
Table 5 shows the low order of accuracy of the present method for the 
Example 3 in comparison with the numerical results presented in Table 
1 and Table 3 for the Example 1 and Example 2, respectively; in which 
the sufficient compatibility conditions are satished. From this one can 
infer that, in practice some of the theoretical compatibility conditions 
seems to be very necessary for high order convergence of the present 
method. Clearly the numerical results presented in Table 1 and Table 
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Table 4. Spectral radius of the transition operator 
for the Example 2. 


e = 2-^ N = U iV = 128 iV = 256 N = 512 N = 1024 
At = 0.5 At = 0.5/4 At = 0.5/4^ At = 0.5/4^ At = 0.5/4^ 


A;=4 

0.01380 

0.68153 

0.90960 

0.97661 

0.99410 

8 

0.25983 

0.72130 

0.93118 

0.98234 

0.99556 

12 

0.32590 

0.77295 

0.93796 

0.98412 

0.99601 

16 

0.34288 

0.78200 

0.94065 

0.98482 

0.99618 

20 

0.35146 

0.78561 

0.94172 

0.98510 

0.99625 

24 

0.35448 

0.78705 

0.94214 

0.98526 

0.99628 

28 

0.35577 

0.78762 

0.94231 

0.98528 

0.99628 

32 

0.35625 

0.78784 

0.94231 

0.98528 

0.99628 


3 verify our theoretical results. 

Previously, the Crank-Nicolson method has been used in the frame¬ 
work of scalar singulary perturbed problem, for instance, in [2] to solve 
one dimensional parabolic problems of convection diffusion type. Re¬ 
cently, Clavero et ah |1] considered the Crank-Nicolson method on 
uniform mesh in time discretization and the central difference scheme 
on standard Shishkin mesh in spatial discretization for a system of two 
coupled time dependent singularly perturbed reaction-diffusion prob¬ 
lems. In this article, to obtain a high order robust approximation we 
considered the Crank-Nicolson method in time direction and a hybrid 
scheme which is a suitable combination of fourth order compact dif¬ 
ference scheme (or HODIE scheme ) and standard central difference 
scheme on a generalized Shishkin mesh in spatial direction. Here it 
is interesting to see how the HODIE technique permits to obtain a 
uniformly convergent method having order bigger than two in spatial 
direction. Earlier, the HODIE scheme for scalar singularly perturbed 
reaction-diffusion problems has been considered in Clavero and Gracia 
[3], and it is proved that the scheme is third order uniformly convergent 
on standard Shishkin mesh. But the extension of new HODIE scheme 
on standard Shishkin mesh is not possible in the case of system of cou¬ 
pled reaction-diffusion problems. It can be seen that the coefficients 
gf’s in (24)-(25) is not always positive at the transition points, due 
to the fact that standard Shishkin mesh is very anisotropic in nature. 
This shows that the operator in (24)-(25) is not of positive type on 
standard Shishkin mesh. At the moment, when N~^ < sfe we can not 
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Table 5. Maximum error and numerical rate of conver¬ 
gence of the present method with uniform step size At in 
time direction and generalized Shishkin mesh S{L) with 
L = L* in spatial direction for the Example 3. 


£ = 2-^ 

N = 64 iV = 128 N = 256 
At = 0.5 At = 0.5/4 At = 0.5/4^ 

iV = 512 N = 1024 
At = 0 . 5/43 ^ q_5/44 

k=A 

3.32E-02 

7.75E-03 

1.91E-03 

4.76E-04 

9.52E-05 


2.10 

2.02 

2.00 

2.32 


8 

3.29E-02 

7.75E-03 

1.91E-03 

4.76E-04 

9.53E-05 


2.08 

2.02 

2.00 

2.32 


12 

1.67E-02 

2.41E-03 

5.17E-04 

1.24E-04 

2.49E-05 


2.80 

2.22 

2.06 

2.32 


16 

1.67E-02 

2.16E-03 

4.41E-04 

1.24E-04 

2.47E-05 


2.95 

2.29 

1.84 

2.32 


20 

1.67E-02 

2.16E-03 

4.41E-04 

1.24E-04 

2.47E-05 


2.95 

2.29 

1.84 

2.32 


24 

1.67E-02 

2.16E-03 

4.41E-04 

1.24E-04 

2.47E-05 


2.95 

2.29 

1.84 

2.32 


28 

1.67E-02 

2.16E-03 

4.41E-04 

1.24E-04 

2.47E-05 


2.95 

2.29 

1.84 

2.32 


32 

1.67E-02 

2.16E-03 

4.41E-04 

1.24E-04 

2.47E-05 


2.95 

2.29 

1.84 

2.32 


EN,At 

3.32E-02 

7.75E-03 

1.91E-03 

4.76E-04 

9.52E-05 

pN 

2.10 

2.02 

2.00 

2.32 



hnd a difference scheme of positive type which is high order uniformly 
convergent on standard Shishkin mesh for system of coupled reaction- 
diffusion problems. To avoid this, one can use the central difference 
scheme in the regular region [r, 1 — t] and the fourth order compact 
difference scheme in (0, r) U (1 — r, 1). But this combination gives 
only second order uniformly convergent result. In order to increase 
the order of convergence and to maintain the positivity of the present 
discrete operator in (24)-(25), we consider a generalized Shishkin mesh 
instead of standard Shishkin mesh. The Lemma 4.1 shows that the dis¬ 
crete operator in (24)-(25) on a generalized Shishkin mesh is of positive 
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Table 6. Spectral radius of the transition operator 
for the Example 3. 


e = 2-^ N = U iV = 128 iV = 256 N = 512 N = 1024 
At = 0.5 At = 0.5/4 At = 0.5/4^ At = 0.5/4^ At = 0.5/4^ 


k=4: 

0.42429 

0.88209 

0.95077 

0.98745 

0.99685 

8 

0.58776 

0.88234 

0.96806 

0.99192 

0.99792 

12 

0.59923 

0.88235 

0.96916 

0.99220 

0.99793 

16 

0.59995 

0.88235 

0.96923 

0.99222 

0.99794 

20 

0.60000 

0.88235 

0.96923 

0.99222 

0.99794 

24 

0.60000 

0.88235 

0.96923 

0.99222 

0.99794 

28 

0.60000 

0.88235 

0.96923 

0.99222 

0.99794 

32 

0.60000 

0.88235 

0.96923 

0.99222 

0.99794 


type and the analysis in Section 4 shows that the scheme (24)-(25) is 
almost fourth order uniformly convergent with respect to the pertur¬ 
bation parameter £. Here we also want to point out one more beneht 
of generalized Shishkin mesh over standard Shishkin mesh in the nu¬ 
merical methods presented in [1] and [5] for parabolic reaction diffusion 
systems. It is proved that the numerical methods presented in [1] and 
[5] have almost second order uniform convergence under the theoretical 
relation N~^ ^ CAt, where 0 < g < 1. Note that the theoretical rela¬ 
tion appeared in the analysis when the barrier function technique was 
used to prove the second order convergence of the regular component 
on standard Shishkin mesh. While if we use generalized Shishkin mesh 
instead of standard Shishkin mesh in |1] and [5] then we can claim al¬ 
most second order uniform convergence in spatial variable without any 
theoretical relation by using the same analysis technique. 

7. Conclusions 

We presented a high order parameter-robust numerical method for 
a system of (M > 2) coupled singularly perturbed parabolic reaction- 
diffusion problem (l)-(3). The problem is discretized using the Crank- 
Nicolson method on an uniform mesh in time direction and a suitable 
combination of the fourth order compact difference scheme and the 
central difference scheme on a generalized Shishkin mesh in spatial 
direction. The essential idea in this method is to use a generalized 
Shishkin mesh in order to attain a high order parameter-robust con¬ 
vergence in spatial variable. The hne parts of standard Shishkin mesh 
and generalized Shishkin mesh are identical, but the coarse part of 
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generalized Shishkin mesh is a smooth continuation of the hue mesh 
and is no longer equidistant. Using this fact we proved that the present 
method is second order uniformly convergent in time and almost fourth 
order uniformly convergent in spatial variable, if the discretization pa¬ 
rameters satisfy a non-restrictive relation. Numerical experiments are 
presented to validate the theoretical results and also the results of the 
experiments indicate that the relation between the discretization pa¬ 
rameters is not necessary in practice. 
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